clear all
set more off

// RDD catchment by political alignment

global input "...\replication\data"
global output "...\replication\outputs"

use "$input\com_data.dta", clear

****************** Costmetics ******************
** Reomove Corsica
ren com_res com
drop if substr(com,2,1) == "A"
drop if substr(com,2,1) == "B"
destring com, replace

****************** Get dist border ******************
merge m:1 com using "$input\com_arcgis.dta"
drop if _merge != 3
drop _merge
** Drop those within 30km of coast or national border
keep if dist_nat_border < 0 

****************** Get neighbor region for all com ******************
tostring com, replace
replace com = "0"+com if strlen(com) == 4
merge m:1 com using "$input\loc_i_raw.dta" 
drop if _merge == 2
drop _merge
destring com, replace

****************** Get dist instrumental border ******************
merge m:1 com using "$input\com_dist_instr.dta"
drop if _merge != 3
drop _merge

destring dep_res dep_work, replace
geonear com x_com y_com using "$input\loc_j_raw.dta", neighbors(dep x_dep y_dep) nearcount(2)
gen other_dep = nid1 if nid1 != dep_res
replace other_dep = nid2 if nid1 == dep_res

keep com ind dist_dep dep_res same other_dep dep_work

****************** Get main vars******************
** Labor markets at home and on the other side
egen lab_mkts = sum(ind), by(dep_work)
preserve
collapse (mean) lab_mkts, by(dep_work)
ren dep_work dep
save lab_mkts.dta, replace
restore

* Get size of home lab market
ren dep_res dep
merge m:1 dep using lab_mkts.dta
drop _merge
ren dep dep_res
ren lab_mkts home_lab_mkt

* Get size of lab market behinf the border
ren other_dep dep
merge m:1 dep using lab_mkts.dta
drop _merge
ren dep other_dep
ren lab_mkts other_lab_mkt


********************************************************
****************** Bunching analysis ******************

local bw 50

****************** Gen treatment ******************

replace dist_dep = -dist_dep if same == 0

drop if same == 0 & other_lab_mkt > home_lab_mkt
drop if same == -1 & other_lab_mkt < home_lab_mkt

gen Treat = 0 if other_lab_mkt < home_lab_mkt
replace Treat = 1 if other_lab_mkt > home_lab_mkt

****************** Regressions ******************

tostring dep_res other_dep, replace

ren dep_res dep
merge m:1 dep using  "$input\pol_color_dep.dta"
drop _merge
ren win_left win_left_res
ren dep dep_res

ren other_dep dep
merge m:1 dep using  "$input\pol_color_dep.dta"
drop _merge
ren win_left win_left_work
ren dep other_dep

gen party_sim = abs(win_left_res - win_left_work)

destring dep_res other_dep, replace

keep if abs(dist_dep) <= 25

preserve
keep if party_sim <= 0.4 // below average using commute dataset for threshold

****************** calonico et al ******************
rdbwselect ind dist_dep, c(0) 
local opt_h=e(h_mserd)

qui: rdrobust ind dist_dep ,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_1

qui: rdrobust ind dist_dep ,c(0) h(`opt_h') covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_2

qui: rdrobust ind dist_dep ,c(0) h(2*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
est store ATE_13

qui: rdrobust ind dist_dep ,c(0) h(0.5*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
est store ATE_14

****************** Results in Table ******************
esttab ATE_1 ATE_2, ///
b(%9.3f) se(%9.3f) stats( h_l N_bw, labels("Bandwidth (km)" "Obs. bw." "Adj. R$^2$" AIC Obs. ) ///
fmt( %9.3gc %9.0gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE")
restore

keep if party_sim > 0.4  // above 

****************** calonico et al ******************
rdbwselect ind dist_dep, c(0) 
local opt_h=e(h_mserd)

qui: rdrobust ind dist_dep ,c(0) h(`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_3

qui: rdrobust ind dist_dep ,c(0) h(`opt_h') covs(dep_res other_dep)
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
sum ind if abs(dist_dep) < `opt_h' & dist_dep < 0
sca m_outc = r(mean)
estadd scalar m_outc=m_outc
est store ATE_4

qui: rdrobust ind dist_dep ,c(0) h(2*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
est store ATE_23

qui: rdrobust ind dist_dep ,c(0) h(0.5*`opt_h') 
sca N_bw=e(N_h_l)+e(N_h_r)
estadd scalar N_bw=N_bw
est store ATE_24

****************** Results in Table ******************
esttab ATE_1 ATE_2 ATE_3 ATE_4, tex ///
b(%9.3f) se(%9.3f) stats( h_l N N_bw m_outc, labels("Bandwidth" "Obs" "Obs. bw." "Mean dep. var. below" ) ///
fmt( %9.3gc %9.0gc %9.0gc %9.3gc)) nolines nogaps star(* 0.10 ** 0.05 *** 0.01) ///
keep(RD_Estimate) coef(RD\_Estimate  "ATE") 

